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Abstract. We derive Copernicus's epicycles from New- 
ton's gravitational force law by assuming that a planet's 
orbit is a perturbed circular orbit, with the perturbation 
defined to be co- rotating with the said orbit. We substi- 
tute this orbit expression into Newton's gravitation law 
and showed that the perturbation satisfies the linear part 
of Hill's oscillator equation for lunar motion. We solve 
this oscillator equation using an exponential Fourier se- 
ries and impose the boundary conditions at the aphelion 
and perihelion to derive the Copernicus's formulas for the 
eccentric, deferent, and epicycle. We show that for small 
eccetricity, the Copernican orbit expression also leads to 
Kepler's law of areas for planetary motion. The formal- 
ism we use is the Clifford (geometric) algebra Cl2,o- 

1 Introduction 

In many introductory physics courses, especially those 
dealing with the history and philosophy of science, 
the Copernican model is taught conceptually but not 
mathematically. [1_ And in undergraduate and graduate 
physics courses, it is not even mentioned at all. One 
possible reason is that, unlike in the case of Kepler's el- 
lipse, Copernicus' epicycles has not been rigorously de- 
rived before from first principles, i.e., from Newton's laws 
of motion and gravitation. So our aim in this paper is to 
present this derivation. But before we do so, let us first 
review the Copernican model. 

Copernicus believed that planets orbit around the sun. 
If the orbit of a planet is circular with the sun at the 
center, then the planet's position in complex form is 



= roe"^°\ 



(1) 



where tq and luq are the planet's orbital radius and fre- 
quency, respectively. But because planets sometimes 



move closest to the sun (perihelion) and sometimes far- 
thest (aphelion) , Copernicus displaced the center of the 
planet's circular orbit a little away from the sun to a new 
point called the eccentric, so that the new position f of 
the planet is 

f = r_i+roe'"«*, (2) 

where r_i is the distance of the eccentric from the sun. 
Actually, the eccentric hypothesis in Eq. ^ is not 
Copernicus' original idea but was already known more 
than a thousand years prior by Ptolemy (though he as- 
sumed that the earth is at rest and not the sun as in 
Copernicus). In fact, if we factor out the exponential 
gti^ot jj^ gq^ Q^ TffQ would arrive at Ptolemy's theorem 
applied by Copernicus in his heliocentric theory: 



roe 



iujot 



= (ro+r_ie-*"'")e' 



-~iuJ0t\ ILOQt 



(3) 



In Ptolemaic terms, r_ie~''^°' is called an epicycle and 
Eq. ([3]) is called the eccentric-epicycle equivalence theo- 
rem. (The actual theorem is stated geometrically. [21 [3] ) 
Notice that the theorem essentially states the equivalence 
of the description of the planet's position in the inertial 
frame (left hand side) and in the rotating frame (quantity 
in parenthesis on the right hand side). 

Yet Eq. ^ is still not consistent with the numerical 
data. Ptolemy resolved this problem by assuming that 
the planet's circular orbit is uniform not with respect to 
the orbit's geometric center but on another point called 
the equant[4l [5]. Though this construction saves the 
appearances, Copernicus claimed that the equant goes 
against the idea of uniform circular motion and for him 
this is " not sufficiently pleasing to the mind" [^ . To rem- 
edy this aesthetic difficulty, Copernicus added on top of 
his original circle in the inertial frame another circle with 
twice the frequency [7j. In complex notation, we write 



r-i 



rQe''^"^ + rie^'"^°^ 



(4) 



where r_i = — 3ri. 



In terms of the orbit's semimajor axis and eccentric- 
ity e, the Copernican expression in Eq. ([4]) becomes 



f = A{-e + e 



iujQt 



1 



^2iujat 



(5) 



for ^,1^—1, 2. That is, 



eie2 = e2ei. 



(11) 
(12) 



as given by Gallavotti[S]. Notice that Eq. ([5]) is different 
from that derived from Kepler's eUiptical orbit for small 
eccentricity: [8] 



f = A{e{l - i) + e"'"' + e(l + £)e"'"'''). (6) 

In this paper, our aim is to show that the Copernican 
expression in Eq. ([5]) is a consequence of Newton's grav- 
itational force law. 

We shall divide the paper into four sections. The first 
section is Introduction. In the second section, we shall 
present a brief tutorial on the Clifford (geometric) al- 
gebra Ch^o for the plane[21 [TOl m [H [T3] , which com- 
bines scalars, vectors and imaginary numbers. We shall 
show how the exponential Fourier series are related to 
eccentrics, deferents, and epicycles. [14j [15] In the third 
section, we shall introduce a pertubation in the planet's 
position in the frame co-rotating with the planet's un- 
perturbed circular orbit and substitute the result to the 
vector form of Newton's law of gravitation. We shall 
show that the perturbation in complex form satisfies the 
linear harmonic oscillator equation 



= s -I- 2£tJos — ;^^o(^ + ^*)^ 



whose scalar and imaginary parts are 



= 
= 



Us + 2cjois 



ZuJqXs, 



(7) 



(8) 
(9) 



These equations are the linear part of Hill's equations 
for lunar motion flHlflT] . We shall solve Eq. ([7]) using ex- 
ponential Fourier series and impose the boundary condi- 
tions at the aphelion and perihelion to derive the Fourier 
coefficients of the Copernican orbit. And the fourth sec- 
tion is Conclusions. 



2 Geometric Algebra 

2.1 Vectors and Complex Numbers 

The Clifford (geometric) algebra C/2,0 is an associative 
algebra generated by two vectors ei and 62 that corre- 
spond to the basis vectors along the x— and y— axis in 
the Cartesian coordinate system. The vectors satisfy the 
orthonormality relation 



The first equation algebraically defines ei and 62 as unit 
vectors by setting their squares to unity; the second equa- 
tion defines the vectors as mutually orthogonal by mak- 
ing their product anticommute. 
Let us define the unit bivector 



I = eie2. 



(13) 



From the orthonormality axiom in Eq. (|10p , it is easy to 
see that i is an imaginary number, 

z^ = 61626162 = -61(6262)61 = -6161 = -1, (14) 

that anticommutes with vectors 61 and 62: 



6ii 

62 1 



62 = -«6i, 

-61 == -i62. 



(15) 
(16) 



Notice that right-multiplying £ to a vector rotates it coun- 
terclockwise by 7r/2. 

In general, a vector a in the two-dimensional space 
spanned by 61 and 62 is given by 



where 



a = 02:61 -|- a„62 = 61a = a 61, 



Q,x 1 ^y'^1 



^x (X^l. 



(17) 



(18) 
(19) 



Equations ([T7]) to P^ relates the vector a to the complex 
number a and its complex conjugate a* . 

If vector b = &a;6i -I- 6j;e2, then the product of vectors 
a and b is 



where 



ab = ab + aAb:= a*6, 



a-b = axbx+ayby, 
aAb = {uxby - aybx)i 



(20) 



(21) 
(22) 



^ly^Li ^^ iiu-) 



(10) 



are the scalar (dot) and imaginary (bivector or planar) 
parts of the product ab = a*b. Notice that the magni- 
tude of the wedge product is that of the cross product 
a X b. (Geometrically, we say that a x b is the vector 
perpendicular to the oriented plane aAb, though tech- 
nically, a X b is not defined in C^2,o — only in 01^^). 



2.2 Circles, Epicycles, and Fourier Series we arrive at 



Because i is an imaginary number, then Euler's theorem 
holds: 

e'^ = cos6l + zsin6', (23) 

where is a real number. If we left-multiply Eq. (|23p by 
ei, we get 

eie'" = eiCos6' + e2sin6', (24) 

where we used Eq. (fT5|l . Equation ([24|l states that eie*^ 
is the vector ei rotated counterclockwise by an angle 9 
(assuming that ei points to the right and 62 points up). 



r = ei(fo +rif/>i +f2V'2)- 




>ei 



Fig. 1. The vector r ~ eire 



i(a;i+0) 



The theorem in Eq. ([2l)l enables us to express the po- 
sition r of a point in uniform circular motion as 

r = eire'('^*+'^) = eir cos{Lut + (/))+ e2rsm{Lot + (j}), (25) 

where r is radius, to is the angular frequency, and cj) is 
the rotational phase angle. Another way to express r is 



r = eifijj, 



where 



20 



r = re 

Aujt 



i, 



(26) 

(27) 
(28) 



are the complex radius and rotor (rotation operator) , re- 
spectively. (See Fig. (1)) 

Let ri and V2 be two rotating vectors: 



ri = eifiV3i=eirie*("i*+'^i', 
r2 = eif2V32=eir2e*('^^*+*^). 

If we displace their sum by a vector rg , 

ro = eifo =eiroe**". 



(29) 
(30) 



(32) 



One way to simplify Eq. (|32p is to set lo\ — lo and L02 = 
luj. So using the definition of the rotor V' in Eq. (|28p. we 
get 

r = ei(fo+riV' + f2V'^). (33) 

The zeroth harmonic is the eccentric; the first, the def- 
erent; and the second, the epicycle. In general, we may 
express the position r in time t as an infinite Fourier 
series: 



r = ei 



00 

E 

A; — — ex 



Tk^^ = y^ eire'"('='^*+'^''). (34) 



k — — ex 



Equation ([34]) represents the Copernican ideal of decom- 
posing an orbit as a sum of epicycles with harmonic fre- 
quencies. 

3 Copernican Dynamics 

3.1 Uniform Circular Orbit 

In Newton's law of gravitation, the equation of motion 
of a planet of mass m revolving around the sun of mass 
M is 



r = -GA/- 



(35) 



where r is the position of the planet with respect to the 
the sun at the origin. 

One to solution to Eq. ((35|) is a circular orbit: 



r = ro = eifoV'o = eiroe' 



J(uJo* + 0o) 



(36) 



where rg is the orbital radius, loq is the orbital frequency, 
and 00 is the orbital phase angle. 

To verify this claim, we first take the derivatives in 
time of the position vector r: 



r = cJoeiifo'(/'o = -«woro, 
i^ = -ijjlva. 



(37) 
(38) 



Next, we take the square of the position r by using the 
conjugation theorem in Eq. (|17p : 



r^ = eifo^oeiro^o = ^I^q ^f^^o = flfa = Tq, (39) 

so that |r| — vq as we expect. And finally, we substitute 
Eqs. dST]) to (glD back to Eq. ^ to arrive at 



GM 



(40) 



(31) which is the circular orbit condition. 



3.2 Linear Perturbation Theory 



Let us assume that the solution to Eq. (pS]) may be ex- 
pressed as a sum. of a cucular orbital position Vq and its 
small correction ri: 



r == ro + Ari, 



(41) 



where A is a perturbation parameter that will be set to 
unity later. If we also assume that the perturbation ri 
lies in the same orbital plane as the original circular orbit 
Tq in Eq. p6p and co-rotating with it, then we may write 
ri as 

ri = eis^po, (42) 



where s is a complex function. Hence, 

r = ei(fo +As)-0o. 



(43) 



Taking the first and second time derivatives of the po- 
sition r in Eq. (^5]) . we get 

r = ei{X's + iujo{fo + Xs))4'o, (44) 

r = ei{Xs + 2XiiL>QS — iLjl{'ro + Xsjjtpo. (45) 



4 Copernican Analysis 

4.1 Epicyclical Fourier Series 

We assume that the solution to the orbital harmonic os- 
cillator equation in Eq. (|50p is an exponential Fourier 
series with luq as the fundamental angular frequency: 

oo 

s= Y. ^fcV-o- (51) 

k — — oo 

The time derivatives of s are 

s = zuJo ^ kak%, (52) 

k — — oo 

oo 

§ ^ -Lol Y^ Pttki^^, (53) 

k — — oo 

while the conjugate of s is 

oo oo 

^*= E «^^o''- E ^-k^o- (54) 



fc=-c 



Equation (j45|l provides the expansion of the left side of 
Newton's gravitation law in Eq. ([35|) . 

On the other hand, to rewrite the right side of the 
gravitation law, we need first to take the square of the 
position vector r in Eq. (|^5)l and retain only the terms 
up to first order in A: 

r^ = (ro + As)*(ro + As) « r^ + A(f*s + fos*). (46) 

Raising both sides of Eq. (|46p to —3/2 power and em- 
ploying the binomial theorem, we get 



Substituting Eqs. ([HI]) to ^5^ back to Eq. ^SU^, we 



^,V~^^o^^*°'''^''*^ 



(47) 



where 



flo = e'*°. (48) 

Multiplying Eq. (gT]) by the position r in Eq. (|33]) yields 

r 1 



rro + ei 



^ ^s-^(s + f,ls*)]i,^, (49) 



2ro 



where we retained only the terms up to first order in A. 
Now, substituting Eqs. (P5|) and ([JT]) back to the grav- 
itation law in Eq. (|35p. we arrive at 



(50) 



= s -I- 2hjQS — -Wq(s + fjls*). 



If we set (/)o = (this means that orbit is not tilted, as 
we shall show later), so that r)o = e**^" = 1, we get Hill's 
oscillator equation in Eq. ([7]). 



get 



0- E {{k^ + 'ik+-)au + -f,lal)i^l (55) 

k — — oo 

after factoring out —loq and rearranging the terms. Be- 
cause the rotors -00 ^-re orthonormal in the Fourier sense, 
then Eq. ((55|) holds only if the coefficient of V'o is zero 
for all k: 

Q^{k^ + 2k + ^)ak + h,lal. (56) 

Solving for the coefficient Ofc in Eq. ([56|) . we get 

a.u = ~\vl{k^ + '^k+^-)al. (57) 

Replacing the index k by — /c, 

ak = -\vl{k^-2k+^)a*_^, (58) 

and substituting the result back in Eq. ([ST]) , we arrive at 

= (fc2 + 2fc + ^){k' - 2/c + ^) - ^ = e{e - 1), (59) 

after factoring out a^ and rearranging the terms. Hence, 
/(: = {-1,0,1}. (60) 



Because the values of the index k are hmited by 
Eq. (|60p . then the Fourier series for the perturbation s 
in Eq. ([5T|) simphfies to 



s = a-iipQ ^ + ao + ai-00- 



(61) 



The relationship between the coefficients a_i and ai in 
Eq. ((6T|) may be obtained by setting fc = 1 in Eq. ((57| : 



0-1 



-3%2a*. 



Similarly, the condition for ag is 



ao = -»}o«o- 



This is satisfied in three possible ways: 



(62) 



(63) 



(64) 



Substituting the expression for s in Eq. (|6ip back to 
the position vector expression in Eq. (j43p . we get 



(65) 



L(a-i + (^0 + ao)'0o + ai-^o)- 



Let us count the number of unknowns in this equation. 
The coefficient a_i is related to ai by Eq. ([62)l . The 
angular frequency Wq in tAq = e'"°* is related to the radius 
ro of fo = roTjo = roe**^" by Eq. (|^ . The phase angle 0o 
of fo is related to that of uq by Eq. (|64p . Thus, there are 
five unknowns in Eq. ([65|) : ai^:, aij,, ro, 0o, and ao. 

However, the orbit of a planet in the plane is com- 
pletely specified in two ways: (a) given the position ri 
and the velocity vi at a particular time ii or (b) given 
the positions ri and r2 at their respective times ii and ^2- 
In other words, there are two constraint vector equations 
that are equivalent to four scalar equations for the com- 
ponents. These four equations can only determine four 
unknowns and not five, so one of our unknowns is super- 
fluous and this must be ao because ao = is a possibility 
in Eq. §^. Thus, Eq. ^ reduces to 



(66) 



ei(a_i -l-foV'o + aiV'o)- 



Because of the similarity of Eq. (I66p to Eq. ([55)1 , we rec- 
ognize eia_i as the eccentric, eifofAo as the deferent, and 
eiciiipQ as the epicycle in the Copernican model. 

4.2 Boundary Conditions: Aphelion and 
Perihelion 

Suppose that at i = 0, the planet is at its aphelion 
position Ta at a distance Tq from the sun at a coun- 
terclockwise angle 7 from the positive x— axis; while at 
t = t/2 = t:/luq the planet is at its perihelion position Vp 
at a distance Vp from the sun at a similar angle from the 
negative x— axis. Imposing these boundary conditions on 



the position vector r in Eq. (j66p yields two simultaneous 
equations: 



fa = eir^e*''' =ei(a_i +fo +ai) 



(67) 



-eirpc''' = ei(a_i - fo + ai). 



Factoring out ei from Eqs. (|67|l and ([68|l and using the 
expression for a_i in Eq. ([S^ . we get 

rae'^ = -3fi^„al + fo + ai, (69) 

-rpC*^ = -ifjoal ~ fo + ai, (70) 

which are two simultaneous equations for fo and ai. 

Deferent. To solve for fo, we take the difference of 
Eqs. dMl) and ([TDD to get 



1 



fo = rofio = roe'*" = -(r^ + rp)e''^ , 


(71) 


ro = ^{ra+rp), 


(72) 


'/'o = 7- 


(73) 



so that 



Thus, the radius ro of the deferent circle is the length 
of the semimajor axis of the orbit; the phase angle 0o is 
the angle of inclination of the semimajor axis from the 
x— axis along ei. 

Epicycle. To solve for the coefficient di, we add the 
Eqs. dSZI) and dM]) to obtain 



■^{ra - rp)e'''' = -377odi + di. 



(74) 



Because di = a^; -I- iay cannot be readily isolated, we 
separate the real and imaginary parts of Eq. (f74|) to get 



2(^0 -rp) cos 7 



ai2;(— 3cos27 + 1) 

-I- aiy(-3sin27), (75) 



^(ra-rp) sin 7 = aij;(-3sin27) 



-f ai.y(3 cos 27 + 1), (76) 



where we used the relation (/)o = 7 in Eq. (|73p . Solving for 
the components aix and aiy and combining the results, 
we arrive at 

a, = -\{ra-rp)e'\ (77) 

Equation ([TZll states that the radius Oi of the epicycle di 
is one-fourth the difference between the aphelion distance 
Ta and the perihelion distance rp. Note the negative sign. 
Eccentric. After knowing di, we use the coefficient 
relation in Eq. ([62]) to solve for d_i: 



d-i = -{ra -rp)e'^. 



(78) 



Equation ((78|) states that the length a_i of the eccen- 
tric is three-fourth the difference between the aphelion 
distance r^ and the perihelion distance rp. 



4.3 Copernican Orbit 

We now substitute the expressions a— coefficients in 
Eqs. (HZl) and ([78]) and that of fo in Eq. ^ back to 
the expression for the position r in Eq. (pS|) to get 

r = eie'''{-{ra-rp) + -{ra+rp)i^o-^{ra-rp)i^Q)- (79) 

If the semimajor axis' inclination angle j = (f>o = 0, then 
Eq. jlSl) reduces to 

3 1 - 1 

r = ei(^(ra - r^) + -(r<j +rp)V'o " ^(^a - »^p)V'o)- (80) 

Equation ([50)1 is our desired approximation of a planet's 
orbit around the sun using eccentric, deferent, and epicy- 
cle in terms of the planet's aphelion ra and perihelion Vp. 
(See Fig. (2)) 

If we employ the definitions of the Keplerian semima- 
jor axis A and eccentricity e, 



A 



1 



(ra +rp), 



' a -r I p 

then we may rewrite Eq. ([50)1 as 

3 - 1 - 

r^eiAi-e + ipo-^f^^o)^ 

which is Eq. ([5]). Or in Cartesian coordinates, 

3 1 

X — A{-e + cos{uJot) ~ -€cos{2ujot)), 

y = A{sm{ujot) - -esin{2uJot)). 



(81) 
(82) 

(83) 

(84) 
(85) 




Fig. 2. Compass-and-straightedge plotting of 
a Copernican orbit with eccentric distance ai, 
deferent radius tq, and epicycle radius oi. The 
orbit's eccentricity is e = 1/3. 



The time derivative of the planet's position r in 
Eq. §^ is 

V = eiiiLjA{i[jo — e-^o)- (86) 

Left-multiplying this by r, 

3 „ 1 ., .. 

rv = iu;oA''i-e + ^^^^-ei;o^){i;o-eijl) 

= £^oA2[-ie^-i + (l + ie2) 

+ ^eV'o-^e'V^o'], (87) 

and separating the scalar and imaginary parts of the re- 
sult, we arrive at 



r-v = cJo^^(— esin(a;ot) 4- -e^sin(2a;oi)), 

1 3 

rAv = iLJoA^{{l + -e'^)- -e'^cos{2u;ot)). (89) 




Fig. 3. The position of a planet in Copernican 
orbit from t = to i = r/2 at a time interval 
of t/12. The orbit's eccentricity is e = 1/3. 

For nearly circular orbits, the eccentricity e w (e is 
0.0167 for earth and 0.0068 for Venus). So dropping the 
e^ terms in Eqs. ([88]) and (|89p . we arrive at 

r-v = -ujoA^esiniuJot), (90) 

rAv = iujoA'^. (91) 

The first equation means that the position and velocity 
of a planet are perpendicular at the aphelion {t — 0) and 
perihelion (t = r/2 = tt/wq); the second equation implies 
that the oriented area 



A = r A (vSt) = iujoA^St 



(92) 



swept by the radius vector r for a small interval of time St 
is constant, which is Kepler's second law (we can always 
perform an integral to show the validity of the law for 
large time intervals). (See Fig. (3)) 



5 Summary and Conclusions 

In this paper, we derived the Copernican system of epicy- 
cles from Newton's gravitational force law in vector form 
via linear perturbation theory in Clifford (geometric) al- 
gebra C^2,o of the plane. We assumed that the planet's or- 
bit is a perturbed circular orbit, where the perturbation 
is defined as a vector co- rotating with the original orbit. 
We substituted this expression into Newton's gravitation 
law. Using binomial expansion, we showed that this per- 
turbation may be represented by a complex function s 
that satisfies the linearized form of Hill's equation for lu- 
nar motion. This equation is a linear harmonic oscillator 
with imaginary damping term and an extra forcing term 
that is proportional to the conjugate s* . 

We solved this oscillator equation using exponential 
Fourier series with the frequency loq of the unperturbed 
circular orbit as the fundamental frequency. We showed 
that only three harmonics are allowed: -1, 0, and 1. 
This result makes the planet's position as an expo- 
nential Fourier series with three harmonics: 0, 1, 2 — 
corresponding to the planet's eccentric, deferent, and 
epicycle. We determined the values of the Fourier coef- 
ficients by imposing that the planet is at its aphelion at 
i = and at its perihelion at t = t/2 = tt/wq- And from 
this we derived Gallavotti's expression for the Copernican 
orbit in terms of its semimajor axis A and eccentricity e. 

We also computed the dot and wedge products of 
the planet's position and velocity. We showed that for 
small eccentricity e, the dot product is proportional to 
~ ain(ujot)', the wedge product is constant, itooA^, which 
implies that the planet's position vector sweeps out equal 
areas in equal times, as given by Kepler's second law. 
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